\[G = \prod_{i = 1}^k g_i\]
\[g_i = \frac{|4X_i - 2| + a_i}{1 + a_i}\]
\[V_i = \frac{1}{3(1 + a_i)^2}; \]
\[ S_i = \frac{V_i}{\prod_{i=1}^k \left( V_i + 1 \right) - 1} \]
\[S_{Ti} = V_j\cdot\prod_{i\neq j}^k(1 + V_i)\]
The function for calculating the variances and hence the sensitivity indices are defined below:
options(scipen = 999)
set.seed(100)
# Define the function to calculate V_i
calculate_Vi <- function(a_i) {
return(1 / (3 * (1 + a_i)^2))
}
# Define the function to calculate S_i
calculate_Si <- function(Vi, V_product) {
return(Vi / (V_product - 1))
}
# Define the function to calculate V_Ti
calculate_VTi <- function(Vi, V) {
return(Vi * V)
}
# Define the function to calculate S_Ti
calculate_STi <- function(VTi, V_product) {
return(VTi / (V_product - 1))
}
# Main script
a <- c(0, 2, 10, 30, 99, 99) # Example values for a
k <- length(a)
# Calculate V_i for each a_i
Vi <- sapply(a, calculate_Vi)
# Calculate the product of (1 + V_i) for total variance
V_product <- prod(1 + Vi)
# Calculate V_Ti for each a_i
VTi <- sapply(1:k, function(i) calculate_VTi(Vi[i], prod(1 + Vi[-i])))
# Calculate first-order sensitivity indices S_i
Si <- Si_analytical <- sapply(Vi, calculate_Si, V_product = V_product)
# Calculate total-order sensitivity indices S_Ti
STi <- STi_analytical <- sapply(1:k, function(i) calculate_STi(VTi[i], V_product))
The arbitrary coefficients are selected. The higher the value, the less significant the input parameter’s impact on the output.
a <- c(0, 2, 10, 30, 99, 99) # Example values for a
# a <- c(6.42,6.42,6.42,6.42,6.42,6.42)
k <- length(a)
The results are then printed:
# Output the results
cat("First-order sensitivity indices (Si):\n")
print(round(Si, 5))
cat("\nTotal-order sensitivity indices (STi):\n")
print(round(STi, 5))
The sum must be equal to 1.
sensobol
package.set.seed(2)
k <- length(a)
params <- paste("x", 1:k, sep = "")
R <- 10^3
type <- "norm"
conf <- 0.95
sobol_Fun <- function(X, vector) {
y <- 1
for (j in 1:length(a)) {
y <- y * (abs(4 * X[, j] - 2) + vector[j])/(1 + vector[j])
}
return(y)
}
sensobol are computed for each
increasing power of 2 with lapply function:N_values <- 2^(7:15)
results <- lapply(N_values, function(N) {
mat <- sobol_matrices(N = N, params = params, type = "R")
y <- sobol_Fun(mat, a)
ind <- sobol_indices(Y = y, N = N, params = params, boot = TRUE, type = "norm", conf = 0.95, R = R)
list(N = N,
Si = ind$results %>% filter(sensitivity == "Si") %>% pull(original),
STi = ind$results %>% filter(sensitivity == "Ti") %>% pull(original))
})
results object is converted to two usable
data.frame objects, for an easier use:The convergence plots for each input sensitivity can be plotted:
# Transform the data to be properly plotted
plot_data <- final_results_Si %>%
pivot_longer(cols = -c(coefficients, Si_analytical),
names_to = "N",
values_to = "Si_value",
names_prefix = "Si_") %>%
mutate(N = as.numeric(N))
# Create the plot
Si_ggplot <- ggplot(plot_data, aes(x = N, y = Si_value, color = factor(coefficients))) +
geom_line() +
geom_point() +
geom_hline(aes(yintercept = Si_analytical, color = factor(coefficients)),
linetype = "dashed") +
scale_x_log10(breaks = unique(plot_data$N)) +
labs(title = "Convergence of Si values to analytical values (dashed lines)",
x = "Sample size (N)",
y = "",
color = "Coefficients:") +
theme_minimal() +
theme(legend.position = "right")
ggplotly(Si_ggplot)
# Transform the data to be properly plotted
plot_data <- final_results_STi %>%
pivot_longer(cols = -c(coefficients, STi_analytical),
names_to = "N",
values_to = "STi_value",
names_prefix = "STi_") %>%
mutate(N = as.numeric(N))
# Create the plot
STi_ggplot <- ggplot(plot_data, aes(x = N, y = STi_value, color = factor(coefficients))) +
geom_line() +
geom_point() +
geom_hline(aes(yintercept = STi_analytical, color = factor(coefficients)),
linetype = "dashed") +
scale_x_log10(breaks = unique(plot_data$N)) +
labs(title = "Convergence of STi values to analytical values (dashed lines)",
x = "Sample size (N)",
y = "",
color = "Coefficients:") +
theme_minimal() +
theme(legend.position = "right")
ggplotly(STi_ggplot)
Pass the cursor over the dots the read the sensitivity values.
Given the simplicity of the function, the convergence is reached quickly, even with a small sample size. These plots can thus be readily used for more complex functions, where convergence might be slower and more interesting to assess.